home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / rng / tt.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  3.8 KB  |  144 lines
  1. /* rng/tt.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_rng.h>
  23.  
  24. /* This is the TT800 twisted GSFR generator for 32 bit integers. It
  25.    has been superceded by MT19937 (mt.c). The period is 2^800.
  26.  
  27.    This implementation is based on tt800.c, July 8th 1996 version by
  28.    M. Matsumoto, email: matumoto@math.keio.ac.jp
  29.  
  30.    From: Makoto Matsumoto and Yoshiharu Kurita, "Twisted GFSR
  31.    Generators II", ACM Transactions on Modelling and Computer
  32.    Simulation, Vol. 4, No. 3, 1994, pages 254-266. */
  33.  
  34. static inline unsigned long int tt_get (void *vstate);
  35. static double tt_get_double (void *vstate);
  36. static void tt_set (void *state, unsigned long int s);
  37.  
  38. #define N 25
  39. #define M 7
  40.  
  41. typedef struct
  42.   {
  43.     int n;
  44.     unsigned long int x[N];
  45.   }
  46. tt_state_t;
  47.  
  48. static inline unsigned long int
  49. tt_get (void *vstate)
  50. {
  51.   tt_state_t *state = (tt_state_t *) vstate;
  52.  
  53.   /* this is the magic vector, a */
  54.  
  55.   const unsigned long mag01[2] =
  56.   {0x00000000, 0x8ebfd028UL};
  57.   unsigned long int y;
  58.   unsigned long int *const x = state->x;
  59.   int n = state->n;
  60.  
  61.   if (n >= N)
  62.     {
  63.       int i;
  64.       for (i = 0; i < N - M; i++)
  65.     {
  66.       x[i] = x[i + M] ^ (x[i] >> 1) ^ mag01[x[i] % 2];
  67.     }
  68.       for (; i < N; i++)
  69.     {
  70.       x[i] = x[i + (M - N)] ^ (x[i] >> 1) ^ mag01[x[i] % 2];
  71.     };
  72.       n = 0;
  73.     }
  74.  
  75.   y = x[n];
  76.   y ^= (y << 7) & 0x2b5b2500UL;        /* s and b, magic vectors */
  77.   y ^= (y << 15) & 0xdb8b0000UL;    /* t and c, magic vectors */
  78.   y &= 0xffffffffUL;    /* you may delete this line if word size = 32 */
  79.  
  80.   /* The following line was added by Makoto Matsumoto in the 1996
  81.      version to improve lower bit's correlation.  Delete this line
  82.      to use the code published in 1994.  */
  83.  
  84.   y ^= (y >> 16);    /* added to the 1994 version */
  85.  
  86.   state->n = n + 1;
  87.  
  88.   return y;
  89. }
  90.  
  91. static double
  92. tt_get_double (void * vstate)
  93. {
  94.   return tt_get (vstate) / 4294967296.0 ;
  95. }
  96.  
  97. static void
  98. tt_set (void *vstate, unsigned long int s)
  99. {
  100.   tt_state_t *state = (tt_state_t *) vstate;
  101.  
  102.   const tt_state_t init_state =
  103.   {0,
  104.    {0x95f24dabUL, 0x0b685215UL, 0xe76ccae7UL,
  105.     0xaf3ec239UL, 0x715fad23UL, 0x24a590adUL,
  106.     0x69e4b5efUL, 0xbf456141UL, 0x96bc1b7bUL,
  107.     0xa7bdf825UL, 0xc1de75b7UL, 0x8858a9c9UL,
  108.     0x2da87693UL, 0xb657f9ddUL, 0xffdc8a9fUL,
  109.     0x8121da71UL, 0x8b823ecbUL, 0x885d05f5UL,
  110.     0x4e20cd47UL, 0x5a9ad5d9UL, 0x512c0c03UL,
  111.     0xea857ccdUL, 0x4cc1d30fUL, 0x8891a8a1UL,
  112.     0xa6b7aadbUL}};
  113.  
  114.  
  115.   if (s == 0)    /* default seed is given explicitly in the original code */
  116.     {
  117.       *state = init_state;
  118.     }
  119.   else
  120.     {
  121.       int i;
  122.  
  123.       state->n = 0;
  124.  
  125.       state->x[0] = s & 0xffffffffUL;
  126.  
  127.       for (i = 1; i < N; i++)
  128.     state->x[i] = (69069 * state->x[i - 1]) & 0xffffffffUL;
  129.     }
  130.  
  131.   return;
  132. }
  133.  
  134. static const gsl_rng_type tt_type =
  135. {"tt800",            /* name */
  136.  0xffffffffUL,            /* RAND_MAX */
  137.  0,                    /* RAND_MIN */
  138.  sizeof (tt_state_t),
  139.  &tt_set,
  140.  &tt_get,
  141.  &tt_get_double};
  142.  
  143. const gsl_rng_type *gsl_rng_tt800 = &tt_type;
  144.